Interplay of shear and bulk viscosity in generating flow in heavy-ion collisions 
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We perform viscous hydrodynamic calculations in 2+1 dimensions to investigate the influence of 
bulk viscosity on the viscous suppression of elliptic flow in non-central heavy-ion collisions at RHIC 
energies. Bulk and shear viscous effects on the evolution of radial and elliptic flow are studied with 
different model assumptions for the transport coefficients. We find that the temperature dependence 
of the relaxation time for the bulk viscous pressure, especially its critical slowing down near the 
quark-hadron phase transition at T c , partially offsets effects from the strong growth of the bulk 
viscosity itself near T c , and that even small values of the specific shear viscosity rj/s of the fireball 
matter can be extracted without large uncertainties from poorly controlled bulk viscous effects. 
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I. INTRODUCTION 

A question of widespread interest is that of the spe- 
cific shear viscosity (shear viscosity per entropy density 
T}/s) of the quark-gluon plasma (QGP) created in nu- 
clear collisions at the Relativistic Heavy Ion Collider 
(RHIC). Ideal (i.e. inviscid) fluid dynamics has been 
quite successful in describing the transverse momentum 
(pr) spectra and the elliptic flow coefficient i>2(pt) of 
the bulk of the thousands of particles created in each, 
say, Au+Au collision [1] . The agreement between theory 
and experiments improves further when one interfaces 
a (3+l)-dimensional ideal fluid description of the QGP 
phase with a hadron cascade during the late expansion 
stage, in order to properly account for the highly viscous 
evolution after hadronization of the QGP [3]. This suc- 
cess strongly suggests that the QGP fireball created at 
RHIC thermalizes very quickly and behaves like an al- 
most perfect liquid [4], which implies that it must be a 
strongly coupled plasma [5-7]. 

On the other hand, the quantum mechanical uncer- 
tainty relation places a fundamental lower bound on 
the specific shear viscosity of any medium [8], and ex- 
plicit computation in a large class of very strongly cou- 
pled quantum field theories (unfortunately not including 
QCD) suggests that this limit is close to the so-called KSS 
bound f | Kgs = 47 ~ 0.08 [9, 10]. While this is a very 
small number (almost two orders of magnitude smaller 
than that of any other known (real) fluid [10, 11], with the 
possible exception of strongly interacting systems of ul- 
tracold fermionic atoms near the unitarity limit [12, 13]), 
it is known that the anisotropic elliptic flow generated in 
non-central relativistic heavy-ion collisions is very sensi- 
tive to shear viscosity [14, 15]. The roots of this sensitiv- 
ity lie in the exceedingly rapid expansion of the heavy- 
ion collision fireballs, especially during the early expan- 
sion stage which is characterized by very large compo- 
nents of the velocity shear tensor [8]. Recent progress 



in performing causal relativistic hydrodynamical simula- 
tions of viscous fluids in 2+1 dimensions [16-28] revealed 
that even very small specific shear viscosities, close to the 
KSS bound, should leave easily identifiable experimental 
signatures, in particular through a suppression of ellip- 
tic flow. That in the experimental data this suppression 
in not large enough to lead to immediate failure of the 
ideal fluid approach suggests that the QGP viscosity, in 
the temperature region probed by RHIC collisions, must 
also be close to the KSS bound [29-33]. 

Viscous hydrodynamics, in comparison with experi- 
mental data, allows in principle for an accurate deter- 
mination of r\j s. In practice, this requires excellent con- 
trol of several other inputs that either are presently not 
known with sufficient accuracy or have not yet been cor- 
rectly implemented in the numerical simulations [34]. 
The largest prevailing uncertainty is related to the initial 
source deformation that drives the elliptic flow which is 
presently not known to better than 20-30% [2, 22, 35, 36] 
(see, however, recent suggestions to eliminate this error 
source [37, 38]). As shown in [22], this leads at present to 
an O(100%) uncertainty in the extracted r\j s value. Two 
other effects of similar magnitude which, however, work 
against each other and may largely cancel, are strong 
viscous effects [2] and the non-equilibrium chemical com- 
position [39-43] in the late hadronic phase. Finally (and 
this is the point of the present paper) bulk viscous ef- 
fects must be included in any study that aims to extract 
the specific shear viscosity [34, 44]. However, even when 
making generous allowances for all these uncertainties, it 
appears clear that the effective shear viscosity to entropy 
ratio of the QGP, averaged over the expansion history 
of the fireballs created in RHIC collisions, cannot exceed 
the conservative upper limit 
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This makes the QGP the most perfect liquid ever ob- 
served in the laboratory. 

In the present paper we use the (2+l)-dimcnsional vis- 
cous relativistic fluid dynamic code VISH2+1 [18, 20, 21] 
to study the effects of bulk viscosity and their interplay 
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with shear viscosity in the buildup of radial and elliptic 
flow. Some preliminary results were reported in [34, 44] 
(see also [45] for related work). Our work is preceded 
by three (0+l)-dimensional studies for systems undergo- 
ing boost-invariant longitudinal expansion without trans- 
verse flow [46-48] which explored the suggestion by Tor- 
rieri, Tomasik and Mishustin [49] that in rapidly expand- 
ing fireballs bulk viscosity can lead to such large negative 
bulk pressures that the fluid becomes mechanically un- 
stable against clustering and cavitation. Since bulk vis- 
cosity is expected to be maximal near the quark-hadron 
phase transition (see discussion in Section 11), those stud- 
ies predicted that bulk viscous effects become important 
mostly during the second half of the fireball expansion 
when the QGP undergoes hadronization. At that time 
the scalar expansion rate 9 = 8^^ (where u^{x) the flow 
4- velocity), which for f -dimensional boost-invariant ex- 
pansion equals 6 = 1/t (where r = ^Jt 2 —z 2 is the longitu- 
dinal proper time, with z indicating the longitudinal or 
beam direction), is already small enough to significantly 
temper the growth (in magnitude) of the (negative) bulk 
pressure, leading to instability problems only for rela- 
tively large peak values for the specific bulk viscosity Q/s 
[46-48]. 

Our work improves on these analysis by including a re- 
alistic initial transverse density profile and the resulting 
transverse flow in the fireball. This has two important 
consequences: (i) The transverse flow increases the ex- 
pansion rate 9, leading to larger bulk pressures | IT | for 
given Q/s. (ii) Some of the matter near the dilute trans- 
verse edge of the fireball experiences large bulk viscosi- 
ties already at very early times where the expansion rate 
9 <~ 1/t is big. This leads to much more severe problems 
with mechanical instability in VISH2+1 than for simple 
f-dimensional expansion, and to correspondingly smaller 
values for the upper limit for (/s that allows for stable 
hydrodynamic evolution. Even more restrictive than the 
condition for mechanical stability is the self-consistency 
constraint for the validity of viscous hydrodynamics it- 
self: the entire framework, which is based on a near- 
equilibrium expansion, breaks down when viscous cor- 
rections to the local equilibrium distribution function be- 
come comparable to the thermal equilibrium terms. We 
will show that this happens, for particles with typical 
momenta p~ 3T, even before the effective total pressure 
becomes negative and mechanical instability sets in [50]. 
While the formalism may be able to qualitatively indicate 
where and when cavitation sets in [46-48] , we doubt that 
the phenomenon itself can be self-consistently described 
within the existing viscous hydrodynamic frameworks. 

Obviously, viscous hydrodynamics can predict the vis- 
cous suppression of elliptic flow reliably only within its 
domain of validity. We therefore restrict our attention 
to the parameter range where the bulk viscous pressure 
stays everywhere sufficiently small that stable hydrody- 
namic evolution is ensured. Within that range (which 
we determine), we study the effects of bulk viscosity and 
of the microscopic relaxation time for the bulk viscous 



pressure on radial and elliptic flow, with and without ad- 
ditional shear viscosity. For a fluid with constant specific 
shear viscosity 2 = j- we find that, depending on initial 
conditions and details of the temperature dependence of 
the relaxation time, bulk viscosity increases the viscous 
suppression of «2(pt) by 5 - 50%. This large range indi- 
cates not only that bulk viscosity is a potentially serious 
contaminant in the extraction of the specific shear vis- 
cosity r//s from elliptic flow data, but also that a robust 
theoretical effort is needed to better constrain the range 
of reasonable values for the bulk viscosity and its associ- 
ated relaxation time. We find that the uncertainty range 
is drastically reduced, to the 10 — 20% level, if we im- 
pose proportionality between the specific bulk viscosity 
and its associated relaxation time, as indicated by kinetic 
theory. The critical growth of the specific bulk viscosity 
near the quark-hadron phase transition is then accom- 
panied by critical slowing down of the dynamics of the 
viscous bulk pressure. This diminishes the bulk viscous 
contribution to the viscous suppression of elliptic flow. 

The paper is organized as follows: In Section II wc 
review the present state of knowledge of the tempera- 
ture dependence of the specific bulk and shear viscosi- 
ties C/s and rj/s and their associated microscopic relax- 
ation times. Based on this analysis we introduce specific 
parametrizations for (£/s)(T) and the bulk pressure re- 
laxation time r n (T) which we use later in the numerical 
simulations. Section III gives a brief summary of specific 
features of the viscous hydrodynamic equations solved in 
this work, referring to earlier work for a more general de- 
scription. In Section IV we discuss generic effects of bulk 
and shear viscosity on the hydrodynamic evolution of fire- 
ball eccentricity and flow and their implications for the 
final hadron spectra and elliptic flow. Section V discusses 
the sensitivity of bulk viscous effects to the initialization 
of the bulk viscous pressure and to its relaxation time. In 
Section VI we explore the range of bulk viscosities that 
allows for stable viscous hydrodynamic evolution. Con- 
sequences of bulk viscous effects for the extraction of the 
specific shear viscosity rj/s from experimental elliptic flow 
data are discussed in Section VII before summarizing our 
findings in Section VIII. 



II. VISCOSITIES AND RELAXATION TIMES 

The present state of knowledge of the viscous proper- 
ties of strongly interacting matter at high temperatures 
is nicely reviewed in [51] to which we refer for details. 
Kinetic theory [52] and experiment [10, 11] show that 
for non-relativistic fluids the specific shear viscosity rj/s 
typically reaches a minimum near the liquid-gas phase 
transition, rising both towards lower temperatures in the 
liquid phase and towards higher temperatures in the gas 
phase. Lattice QCD [53], perturbative QCD [54], and 
hadron cascade simulations [55] indicate that relativis- 
tic QCD matter behaves analogously, but with the liquid 
and gas phases interchanged (the liquid QGP phase ex- 
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ists at higher temperature than the hadronic gas phase) . 
According to perturbative [54] and lattice [53] QCD, the 
increase of 77/s with temperature in the QGP phase is 
weak over the temperature range explored in RHIC col- 
lisions, suggestion the use of a constant 77/s for the QGP 
in hydrodynamic simulations. We here use the smallest 
value for this constant permitted by the KSS conjecture 
[10], f| KSS = ^ w 0.08, in order to extract reasonable 
upper bounds for the uncertainties introduced by bulk 
viscous effects into the extraction of such a small value 
from experimental data. [We note that the assumption of 
a constant 77/s is unacceptable for quantitative attempts 
to extract it from heavy-ion collision data; at least one 
must account for a significant increase of this ratio during 
hadronization and in the late hadronic phase [55] .] 

The relaxation time r w for the shear viscous pressure 
tensor tt^ v has been computed for a relativistic Boltz- 
mann gas [56, 57], in weakly coupled QCD [58], in lat- 
tice QCD [59] , and in Af = 4 SYM theory at infinite cou- 
pling [60—62]. The results can be presented in the form 
r rr = ^^, with A bracketed by 2<A<6. We here use 
t — 1 2 — 1A. 

Ts 2ttT- 

Theoretical knowledge of the specific bulk viscosity (/s 
is more murky. For a non-interacting system of massless 
quanta it vanishes exactly, due to conformal invariance. 
Interactions lead to deviations from zero that usually re- 
main small, except near phase transitions where the sys- 
tem may develop large correlation lengths [63-65]. Ki- 
netic theory gives - = k (|— c^) where k = 5/3 in relax- 
ation time approximation [67] and k = 15 for a system of 
photons radiated by massive particles in thermal equilib- 
rium [68]. The complete leading order result for weakly 
coupled QCD [69] is roughly consistent with the latter 
of these two, but adds a weak (decreasing) temperature 
dependence. At high temperatures c^~^, so the ratio 
£ is small of second order in the deviation. For strongly 
coupled N = 4 SYM theory one obtains a lower bound 
for this ratio which is only linear in this deviation and 
thus much larger: - > 2 (|— c^) [70]. For the hadron gas 

different authors [71-73] agree that ^ -C 1 just below T c 
and that this ratio decreases towards lower temperatures. 
There is no agreement on the sign of the temperature 
dependence of the specific bulk viscosity Q/ s itself which 
according to [71] increases towards lower temperature for 
massive pions, but decreaes for massless pions [73]. How- 
ever, there are general arguments [64, 74] that support 
the idea that (/s should peak near the quark-hadron 
phase transition, due to long-range correlations related 
to the restoration of chiral symmetry; at a second-order 
critical point £/s is predicted to diverge [65, 66]. 

We assume here that (/s quickly approaches zero once 
T decreases below T c ; above T c , we parametrize it as 
§ = 27(5 — c s) (which corresponds to the Buchel bound 
[70] for 77/s = l/(47r)), using c 2 s {T) extracted from the 
same lattice QCD data [75] that we used for our equa- 
tion of state (EOS L, see [21] for details). The factor 
(| — cf) increases as we approach T c from above; the 
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FIG. 1: (Color online) Parametrization of the specific bulk 
viscosity (/sasa function of temperature. (C is a multiplica- 
tive scaling factor for the entire function, see text.) 

resulting increase of Q/ s is qualitatively, but not quanti- 
tatively consistent with a direct extraction of Q/s from 
lattice QCD [76] (for a critical discussion of this extrac- 
tion see [65] ) and with recent work in "holographic QCD" 
[77] . We connect our parametrization above T c to the as- 
sumed zero value for £/s well below T c by interpolating 
with a Gaussian function (see Fig. 1). This results in a 
peak value (£/s)(T c ) ~ 0.04 - about half as big as our 
choice for the (temperature independent) specific shear 
viscosity 77/s and consistent with strong coupling esti- 
mates for strongly coupled conformal field theories using 
the AdS/CFT correspondence [78] (which was also used 
by Buchel when deriving his bound) and with holographic 
QCD [77]. It is, however, more than 10 times smaller 
than both the lattice QCD value extracted by Meyer [76] 
and a recent AdS/CFT-based estimate by Buchel for a 
non-conformal plasma [66]. We will see that this factor 
10 has crucial implications for the applicability of viscous 
hydrodynamics. To simulate larger bulk viscosities, we 
scale the function shown in Fig. 1 (to which we will refer 
as "minimal bulk viscosity" for brevity) by a constant 
factor C > 1. 

Finally, we must specify the relaxation time r n for the 
bulk viscous pressure II about which even less is known 
theoretically. In Israel-Stewart theory (both in its macro- 
scopic and microscopic kinetic formulation [56]) one has 
T n = (f3o where /3o is some combination of thermal equi- 
librium integrals. This suggests that, if £/s peaks near 
T c due to growing correlation lengths, so does the re- 
laxation time r n for the bulk pressure ("critical slowing 
down" [79]). Buchel [66] found that in theories where 
the specific heat diverges at T c , cy <~ l/^/|l— T c /T\, the 
relaxation time can actually diverge at T c even if (/s 
remains finite, i.e. as T— >T C , r n grows more strongly 
than (/s (see also the discussion in [65]). We use the 
parametrization 
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with t = 120fm/c. (1) 
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This increases linearly with £/s as T^T C , but imposes 
a non-zero lower bound on r n , for reasons of numeri- 
cal stability of VISH2+1. For comparison we also study 
two constant relaxation time values, r n =0.5 and 5fm/c, 
roughly corresponding to the smallest and largest values 
of Eq. (1) for temperatures 1 < T/T c < 2 if we set C = 1. 

III. VISCOUS HYDRODYNAMICS 

We solve the following second order viscous hydrody- 
namic equations ( "Israel-Stewart equations" [56, 80-82]), 

£^7^ = 0, T^ v = eu»u v - (p+n)A^ +tt^, (2) 

cn = -±<n + < S) -infM^), <4) 

in the two transverse spatial directions and time ((2+1)- 
d), implementing boost-invariant longitudinal expansion 
along the beam direction. We assume zero net baryon 
density and thus vanishing heat conductivity. Here, T^ v 
is the energy momentum tensor, 7r Ml/ is the shear pressure 
tensor, and II is bulk pressure. cZ M denotes the covariant 
derivative components (see [16, 24] for details) in the 
curvilear coordinates (t, x,y,rj a ) where r = Vt 2 —z 2 is 
the longitudinal proper time and rj s — ^ In is the 
space-time rapidity. A Ml/ = g^ h '—u fJ -u L ' projects onto 
the spatial components in the local rest frame (here 
qiiv _ diag(l, — 1, — 1, — 1/t 2 ) is the metric tensor in 
(t, x, y, 77s) coordinates); \7 tl — A ,ll/ d u is the spatial gra- 
dient and D = u^d^ is the time derivative in that frame. 
The driving forces for the shear and bulk viscous pres- 
sures are the (symmetric and traceless) velocity stress 
tensor o^ v = V^u u > = ±( W+ W) - \A^ U B and the 
scalar expansion rate 9 = d v u v = V v u", respectively. The 
shear and bulk viscosities 77 and £ and their associated re- 
laxation times Ttt and m were discussed in the preceding 
Section. 

The last terms in Eqs. (3) and (4) are of second order 
in deviations from local equilibrium. For conformal sys- 
tems they can be written in various equivalent forms, up 
to higher order corrections [21]. Even for non-conformal 
systems, such as QCD with the equation of state EOS L 
used here (see below), the difference between the terms 
as written down here and their various conformal approx- 
imations are numerically insignificant [21] unless incon- 
sistently large relaxation times are used. Other second 
order terms that should be allowed for on the right hand 
sides of Eqs. (3) and (4) were identified in [60, 83], and 
some of their coefficients were derived in the weak cou- 
pling limit in [58]. Recent code verification efforts by 
the TECHQM Collaboration [84] indicate that these ad- 
ditional terms have very little numerical influence. We 
therefore ignore them in the present study. 



The explicit form of Eqs. (2-4) for longitudinally boost 
invariant (i.e. f^-independent) systems is given in [16]. 
The equations are closed by providing an equation of 
state for which we use EOS L as described in Ref. [21]. 
We study Au+Au collisions with the same initial con- 
ditions for the starting time To = 0.6fm/c and Glauber 
model initial energy density profiles as in Ref. [21], 
with peak density e = e(r=0, r ; 6=0) = 30 GeV/fm 3 in 
central (6=0) collisions. For the viscous pressures 
we use either zero (U(x, y, t ; 6) = n fJ - l, (x, y, t ; 6) = 0) or 
Navier-Stokes initial conditions (II (x, y, t ; b) — — C^( T o) 
and TT ft,/ (x, y, To; b) = 2r\o^ v (to)) , calculated from the ini- 
tial velocity profile (which does not depend on x, y and 6, 
due to the absence of initial transverse flow) . The actual 
choice will be noted when discussing the results. As in 
[21] we end the hydrodynamic evolution and compute the 
final hadron spectra on a freeze-out surface of constant 
temperature Td ec = 130 MeV. 



IV. VISCOUS EVOLUTION AND SPECTRA: 
GENERIC FEATURES 

In this section we compare generic effects on the hy- 
drodynamic evolution and final particle spectra caused 
by shear and bulk viscous effects separately. (Their com- 
bined effects will be explored in Sects. V-VII.) To this 
end we perform hydrodynamic comparison runs for cen- 
tral (6 = 0) and non-central (6 = 7fm) Au+Au collisions, 
using identical initial and freeze-out conditions, for (i) an 
ideal fluid, (ii) a viscous fluid with only minimal shear 
viscosity 7 = 35: , and (iii) a viscous fluid with only "min- 
imal bulk viscosity" (C = 1) as defined in Sec. II. In 
the viscous runs we choose Navier-Stokes initial condi- 
tions for the viscous pressures and equal relaxation times 
Tjr = r n = (As a caveat we note that in the bulk vis- 
cous case results can depend sensitively on the initial 
conditions, depending on the characteristics of the relax- 
ation time for the bulk viscous pressure - see discussion 
in Sec. V.) The results for case (ii) supplement those for 
the smaller Cu+Cu collision system studied in [18, 20] 
(although for a more realistic equation of state and using 
Eq. (3) instead of the "simplified Israel-Stewart equa- 
tion" employed in those earlier papers). 



A. Hydrodynamic evolution 

Figure 2a shows the time evolution of the local tem- 
perature in central Au+Au collisions for the three cases. 
(We plot the temperature at a radius r = 3 fm from the 
fireball center since at r = the curves for cases (i) and 
(iii) arc almost indistinguishable.) Compared with the 
ideal fluid, shear viscosity reduces the work done by lon- 
gitudinal pressure and thus slows down the cooling pro- 
cess during the early stage; during the middle and late 
stages, shear viscosity accelerates the cooling since the 
positive transverse shear pressure leads to stronger radial 
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FIG. 2: (Color online) Time evolution of the local tempera- 
ture (a) and average radial flow (b), from ideal hydrodynam- 
ics (dashed blue), viscous hydrodynamics with only minimal 
shear (solid red) or bulk (solid brown) viscosity. In (b) the 
average radial flow is calculated with the Lorentz contracted 
energy density 7±e in the transverse plane as weight func- 
tion. The inset in (a) shows the late evolution with increased 
resolution. 



flow than for the ideal fluid (see Fig. 2b and Ref. [20] for 
a full discussion). At late times, the shear viscous fireball 
thus cools more rapidly than the ideal fluid [18, 20]. 

Bulk viscosity, on the other hand, reduces the work 
done in all three directions, due to the isotropic negative 
bulk pressure II ~ —£6 driven by the positive expansion 
rate 6 > 0. As a result, radial flow develops less rapidly 
than for the ideal fluid (Fig. 2b), and the bulk viscous 
fluid cools (slightly) more slowly than the ideal one dur- 
ing all expansion stages (Fig. 2a). While the expansion 
rate is largest at very early times, the bulk viscosity then 
is very small throughout the fireball, except for a thin 
region near the transverse edge of the fireball where the 
matter is close to T c ; bulk viscous effects are therefore 
almost negligible until most of the matter cools down to 
T c . At this time the longitudinal expansion rate has sig- 
nificantly decreased, but transverse expansion picks up 



some of the the slack, and we see significant bulk viscous 
effects on radial flow evolution between 5 and lOfm/c. 
Surprisingly, the consequences for the cooling rate are 
significantly smaller than in the shear viscous case: with 
the parameters studied here, the cooling rates for the 
ideal and bulk viscous fluids almost agree. 

We now turn to non-central collisions. To describe 
the fireball deformations in configuration and momen- 
tum space, we use the spatial eccentricity e x = ^ 2 ~^ 2 ^ 
(where ((...)) denotes an energy density weighted av- 
erage over the transverse plane [85]) and the momen- 

irpXX _rpyy\ 

turn anisotropies e„ ~ jrfkvTjrbvi (defined in terms of un- 
weighted averages over the transverse plane of compo- 
nents of the ideal fluid part of the energy-momentum 
tensor, and thus measuring only the collective flow 
anisotropy [20] ) and e' p = ^ Txx+TMy ^ (defined in terms 
of the total energy-momentum tensor that contains the 
viscous pressure components and thus additionally in- 
cludes microscopic momentum anisotropies in the local 
rest frame of the fluid [20]). 

Figure 3a shows the time evolution of the spatial ec- 
centricity e x for non-central Au-I- Au collisions at b = 7 fm. 
Compared with the ideal fluid, bulk viscosity decelerates 
whereas shear viscosity initially accelerates the decrease 
with time of the spatial eccentricity e x . This is a direct 
consequence of the weaker radial flow in the bulk viscous 
case and the stronger radial flow in the shear viscous case. 
It is easy to see that isotropic radial expansion is enough 
to decreases the spatial eccentricity e x [85]; anisotropic 
flow, with larger flow velocities in the reaction plane than 
perpendicular to it, only accelerates the rate with which 
it decreases. At late times, the eccentricity for the shear 
viscous fluid decreases more slowly than for the ideal one, 
since the ideal fluid develops stronger elliptic flow (see 
Fig. 3b and following discussion). In contrast, in the bulk 
viscous case the slower rate of decrease of the eccentricity 
is caused by weaker radial flow. 

As the spatial eccentricity of the fireball decreases, 
its momentum anisotropy increases. This is shown in 
Fig. 3b. The dash-dotted lines for e p , which takes 
only the ideal fluid part Tq V — (e+p)u^u v — pg^ u of the 
energy-momentum tensor into account, show how hy- 
drodynamic forces convert the spatial anisotropy into a 
flow anisotropy. We observe that at early times the flow 
anisotropy e p rides on the developing radial flow: com- 
pared to the ideal fluid, e p develops a little faster in the 
shear viscous fluid but a little more slowly in the bulk 
viscous case, following the evolution of radial flow. 

The difference between e p (dash-dotted lines) and e' p 
(solid lines) stems from the viscous pressure components 
in the energy- momentum tensor. It reflects a contribu- 
tion to the momentum anisotropy that does not arise 
from anisotropic collective flow but from viscous devia- 
tions from isotropy of the microscopic momentum dis- 
tribution / = fcq+Sf in the local fluid rest frame [20] , 
accounted for by the non- ideal terms in T Miy . Figure 3b 
shows that for the shear viscous fluid these viscous correc- 
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FIG. 3: (Color online) Time evolution of spatial eccentricity 
e x (a) and momentum anisotropy e p , e' p (b), from ideal and 
viscous hydrodynamics (see text for details). In (b), the differ- 
ent symbols along the bulk viscous fluid lines indicate central 
(r—0) freeze-out times for different freeze-out temperatures 
as described in the legend. 



tions are large and negative (reflecting local momentum 
anisotropics pointing out of the reaction plane [20], i.e. 
opposite to the flow anisotropy), especially at early times 
when the expansion rate and shear velocity components 
are large. In contrast, the bulk viscous fluid shows 
significant viscous corrections only after about 2.5fm/c, 
lasting until about 8fm/c, which is when in these non- 
central (6 = 7 fm) Au+ Au collisions the bulk of the mat- 
ter passes through the phase transition where Q/s is large 
(c./. the temperature markers on the curves shown in 
Fig. 3b). The bulk viscous pressure contribution to e' p 
is positive, i.e. pointing into the reaction plane, parallel 
to the collective flow anisotropy. (This was also recently 
pointed out by Monnai and Hirano [86].) 

For both shear and bulk viscosity, the sign of the vis- 
cous pressure contributions to e' p obeys a "Lenz rule": 
they act against the radial flow driven effects on the mo- 
mentum anisotropy. At late times these contributions be- 
come small in both cases, in the bulk viscous case driven 
by the rapid disappearance of C below T c (as modeled 



by us), in the shear viscous case by the more gradual 
vanishing of the shear velocity tensor [20]. We see 
that in the bulk viscous case the radial flow effect on e' p 
eventually wins over that from the local deviation from 
equilibrium Sf, whereas the opposite is true for the shear 
viscous fluid. In both cases, the net effect at freeze-out is 
thus a viscous suppression of the momentum anisotropy 
below the ideal fluid limit (see purple stars in Fig. 3b). 
The viscous suppression of e' p resulting from shear vis- 
cosity is 4-5 times stronger that that arising from bulk 
viscosity. 



B. Spectra and elliptic flow 

From the hydrodynamic output at decoupling temper- 
ature Xdec the spectra and their azimuthal anisotropy, 
in particular the elliptic flow coefficient V2(pr), are com- 
puted with a modified Cooper-Frye algorithm that takes 
into account that in viscous hydrodynamics the local 
phase-space distribution f(x,p) on the freeze-out surface 
slightly deviates from thermal equilibrium, f = f C q+Sf, 
due to small but non-zero viscous pressure components 
[15, 20, 87, 88]. Figure 4a shows the pion p^-spectra for 
central Au+Au collisions from ideal and viscous hydro- 
dynamics. Compared to the spectrum from ideal fluid 
dynamics, shear viscosity leads to flatter spectra while 
bulk viscosity generates steeper ones. This is a direct 
reflection of the stronger radial flow caused by the posi- 
tive transverse shear pressure and the weaker radial flow 
resulting from the negative bulk viscous pressure. The 
slightly larger normalization of the viscous spectra is a 
consequence of viscous entropy production which leads 
to larger final multiplicities [21, 87]. 

Figure 4a shows that shear and bulk viscosity act 
against each other in how they affect the slope of the pr~ 
spectra. When both viscosities are included together in 
the viscous calculations, this reduces the amount of read- 
justment needed in the initial conditions when fitting the 
measured transverse momentum spectra with viscous in- 
stead of ideal fluid dynamics [87] . The differential elliptic 
flow V2(pt) for soft pions, on the other hand, is affected 
by both bulk and shear viscosity in the same way: Fig- 
ure 4b shows that the viscous reduction of e' p in Fig. 3b 
translates directly into reduced elliptic flow w 2 of the fi- 
nal hadron spectra. This is true in particular in the low- 
Pt region px < 1 GeV/c (see inset in Fig. 4b). At larger 
Pt, the shear viscous «2 suppression further increases, 
due to a negative contribution from Sf along the freeze- 
out surface. In contrast, bulk viscosity increases V2{pt) 
above lGeV/c because it results in steeper px spectra. 
At pr = 0.5GeV/c (approximately the mean transverse 
momentum for pions) we find that minimal bulk viscos- 
ity suppresses vi by ~ 5% while minimal shear viscosity 
leads to a ~ 20% suppression. 

If this were the complete story, the additional ~ 5% 
bulk viscous vi suppression would lead to a ~ 25% re- 
duction of the value for rj/s that one might extract from 
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FIG. 4: (Color online) pr spectra and elliptic flow V2{pr) for directly emitted pions (i.e. without resonance decay contributions). 



experimental elliptic flow data, by comparing them with 
an ideal fluid dynamical baseline as advertised in [22]. 
This is a non-negligible effect. Since the input for the 
bulk viscosity used in our calculations is fraught with 
large theoretical uncertainties, as discussed in Sec. II, 
this points to a likelihood for correspondingly large un- 
certainties in the empirical extraction of rj/s from ex- 
periment. In the following Section we follow this line of 
thought further, by investigating the additional sensitiv- 
ity of bulk viscous effects to the initial conditions for the 
bulk viscous pressure and to its relaxation time. 



V. SENSITIVITY OF BULK VISCOUS 
DYNAMICS TO INITIAL CONDITIONS AND 
RELAXATION TIMES 

In this Section we now focus entirely on bulk viscosity 
and investigate what happens when we change the initial 
value for the bulk viscous pressure II and its relaxation 
time r n . We study only "minimal bulk viscosity" as de- 
fined in Sec. II (i.e. (7=1), leaving the discussion of 
larger values to the following Section. 

Figure 5 shows, for peripheral Au+Au collisions at 
6 = 7 fm, the differential elliptic flow for pions (a) and the 
time evolution of the average value of the bulk viscous 
pressure (II) (b), for the two initial conditions (zero and 
Navier-Stokes) for II and the three choices of relaxation 
time scales r n discussed in Sec. II. At px = 0.5GeV/c, 
Fig. 5a indicates a bulk viscous i>2 suppression that 
ranges (for C = 1) from ~ 2% to ~ 10%. For the short re- 
laxation time r n = 0.5 fm/c (solid and dotted green lines) 
the suppression is insensitive to the initialization of II, 
yielding about 8% suppression below the ideal fluid value 
at px = 0.5 GeV/c for both zero and Navier-Stokes initial 
values. Figure 5b explains the underlying reason for this 
observation: for this short relaxation time, (II) quickly 
loses all memory of its initial value, relaxing for both ini- 
tial conditions to the same trajectory after ~ 1 — 2 fm/c 
(i.e. after a few relaxation times, similar to what we saw 



earlier [20] for the shear pressure components). This also 
demonstrates that most of the finally observed bulk vis- 
cous V2 suppression is generated during the middle part 
of the expansion when most of the matter cools through 
the phase transition. If it were dominated by large nega- 
tive bulk viscous pressures in the outer shell of the fireball 
at early times, we should see stronger sensitivity to the 
initial value for II. 

This changes completely if we chose a 10 times longer 
relaxation time, r n =5 fm/c (solid and dotted magenta 
lines in Fig. 5). Now the bulk viscous «2 suppression be- 
comes extremely sensitive to the initialization of II. For 
zero initialization, the average bulk pressure (II) always 
remains small, leading to very small (0(2%)) final sup- 
pression effects for uq. For Navier-Stokes intialization, 
(II) is initially very large and negative (due to the large 
initial expansion rate) and, instead of relaxing to smaller 
values as predicted by Navier-Stokes theory (and real- 
ized by the solid green line corresponding to short r n ) , it 
remains larger than the N-S value for about 4 fm/c. As 
seen in Fig. 5a, this leads to a much larger V2 suppression 
of about 10%. 

The "critical slowing down" scenario which uses a tem- 
perature dependent relaxation time that follows the be- 
haviour of (T) is shown by the solid and dash-dotted 

black curves in Fig. 5. In this case the bulk viscous pres- 
sure quickly relaxes to its Navier-Stokes value in the in- 
terior of the fireball where the temperature is high and 
the relaxation time is short; near the edge of the fire- 
ball, however, where the temperature is near T c and the 
relaxation time is long, it remembers its initial value (ei- 
ther zero or the large negative initial N-S value) for a 
long period. With some reflection one convinces oneself 
that this implies that for both zero and N-S initializa- 
tions the magnitude of the average bulk viscous pressure 
(II) remains below the value observed for the short re- 
laxation time. This is seen in Fig. 5b when comparing 
the black and green curves. Correspondingly, the viscous 
i>2 suppression seen in part (a) of the Figure is for both 
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FIG. 5: (Color online) (a) Differential elliptic flow V2(pr) for 
directly emitted pions (without resonance decays) from ideal 
and viscous hydrodynamics, including only minimal (C = I) 
bulk viscosity, (b) Time evolution of the bulk pressure (11) 
averaged over the transverse plane (weighted by the energy 
density) from viscous hydrodynamics. Different curves cor- 
respond to different initializations and relaxation times, as 
indicated (see text for discussion). 



initializations smaller for the "critical slowing down" sce- 
nario than for a short constant relaxation time. When 
comparing the "critical slowing down" scenario with the 
long constant relaxation time, the viscous vi suppression 
is significantly smaller for N-S initialization ((9(7%) vs. 
O(10%)) and about equally small (0(2%)) for zero ini- 
tialization. 

Since these findings contradict at least our own naive 
first expectations, we briefly reiterate the main point: 
taking into account the critical slowing down of the bulk 
viscous pressure dynamics near T c where £/s becomes 
large leads to weaker bulk viscous suppression effects on 
the elliptic flow than seen for both short and long con- 
stant (i.e. T-independent) relaxation times r n . 



VI. LARGER C/s AND THE BREAKDOWN OF 
VISCOUS FLUID DYNAMICS 



As noted in Sec. II, the peak value of £/s in our 
parametrization shown in Fig. 1 is about 10 times smaller 
than some other estimates [66, 76]. When one tries to 
simply multiply the function shown in Fig. 1 by C = 10, 
one finds that (except for special circumstances discussed 
below) the viscous hydrodynamic code crashes. The rea- 
son is that sufficiently large bulk viscosity can lead to 
fireball regions where the effective total isotropic pressure 
p+Tl (thermal + bulk viscous pressure) becomes negative 
and the medium becomes mechanically unstable and will 
tend to break up [46-49]. In fact, since certain com- 
ponents of the shear viscous pressure (in particular its 
longitudinal component ir^) are usually also negative, 
instability can set in even somewhat earlier [47, 48]. In 
numerical simulations this manifests itself through the 
exponential amplification of local numerical errors which 
will eventually stop the code from running. 

We point out that even before the fluid becomes me- 
chanically unstable one has left the region of applica- 
bility of viscous hydrodynamics. The viscous hydrody- 
namic formalism is based on a near-equilibrium expan- 
sion; its validity assumes that the viscous corrections to 
the energy-momentum tensor are small compared with 
the ideal fluid terms. In other words, if the condi- 
tion (|II| + |7r A " y |)/(e+p) <C 1 is violated for any compo- 
nent (//f), the evolution based on equations (2-4) can no 
longer be trusted. Ignoring the shear pressure and set- 
ting e+p — sT^iAp for a QGP, the instability threshold 
p+II = translates into |II|/ (e+p) » i which is not suffi- 
ciently small to trust the continued validity of the equa- 
tions. The following alternate consideration leads to the 
same conclusion: If the fluid can be described by quasi- 
particles, the viscous terms in the energy-momentum ten- 
sor correspond to deviations of the local phase-space dis- 
tribution f(x,p) = feq+Sf from local equilibrium. Using 
Grad's 14-momcnt method, the deviation Sf is expanded 
up to quadratic order in momentum [15, 56, 57, 83, 86] 
and (for a fluid with only bulk viscosity and massless 
particles at midrapidity y = 0) can thus be written in the 
form 



Pt 



n 



si _ 

f C q " T 2 e+p ' 



(5) 



where a is a slowly varying function of temperature with 
magnitude of order unity [86] . When p+TJ = such that 
— — j , this means that for midrapidity particles with 
typical thermal momenta pr — 3T the deviation 5f/f cq 
is negative with magnitude 1 or larger, rendering the to- 
tal distribution function / negative, which is unphysi- 
cal. Clearly the deviations from local equilibrium are 
too large and the formalism breaks down. 

In this Section we explore the range of bulk viscosities 
that are allowed without leaving the region of validity of 
second-order (Israel-Stewart) viscous hydrodynamics. As 
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FIG. 6: (Color online) Upper limits for f/s (a) and viscous 
entropy production (b) as a function of bulk viscous relax- 
ation time T n , for zero and Navier-Stokes initialization. The 
stars indicate the results for the temperature dependent re- 
laxation time (1) with Navier-Stokes initial conditions, for 
To = 0.6 fm/c. For r n given by Eq. (1) and zero initial con- 
ditions, there is no upper limit for £/s, i.e. the fiud remains 
stable for all values of C. 

in the preceding Section, we study both zero and Navier- 
Stokes initial conditions and the same three choices for 
the bulk viscous relaxation time r n , but we now also in- 
clude runs where the fluid has an additional shear viscos- 
ity r)/s= (l-^2)/(47r), with shear viscous relaxation time 
= 3rj/(sT), and we vary the time To when we start the 
hydrodynamic evolution. (For later starting times, we 
downscale the initial peak entropy density sq such that 
the total initial entropy ~ sqTq is held constant.) For 

the specific bulk viscosity f |J (T) we take the functional 

form shown in Fig. 1, but multiplied by an arbitrary con- 
stant C > 1. For each set of initial conditions, r n , and r\j s 
we determine the largest value C max that still allows for 
stable running of the code, i.e. where the effective to- 
tal isotropic pressure p+II does not violate the stability 
criterium p+II > anywhere inside the freeze-out surface. 

Figure 6 shows the upper limit (C/s)max(Ic) (on the 
left vertical axis) and the corresponding maximal C-value 
Cmax (on the right vertical axis) as a function of the 
bulk viscous relaxation time r n . We see that it depends 
strongly on the initialization. 

For Navier-Stokes (N-S) initial conditions ((/s) ma , x (T c ) 
is insensitive to the relaxation time t_ . In this case the 



magnitude of the average bulk pressure II decreases more 
or less monotonically with time (see Fig. 5b). Violations 
of the positivity condition p+Tl > thus always happen 
at the starting time To, at transverse positions where 
the matter is close to the phase transition. This leads 
to a (C/s)max(Tc) that is controlled by initial conditions 
and independent of the relaxation time. (This includes 
the temperature dependent relaxation time (1) - see the 
star in Fig. 6a.) Correspondingly, (C/s) m ax(?c) does not 
depend on the value of 77/s when shear viscosity is in- 
cluded. When one starts the hydrodynamic evolution 
later, (C/s) max (T c ) increases with To. The dependence 
on To arises from the strong dependence of the initial 
bulk pressure II = —£0 = —(,&■ u on To, through the ex- 
pansion rate 0(tq) = 1/tq. This is illustrated by the solid 
red, dashed magenta and dotted orange lines in Fig. 6a: 
as one increases To from 0.6 to 1 and 2fm/c, the maximal 
(C/ s )max(T'c) increases from 0.05 to 0.09 and 0.18. 

For zero initialization II(to) = one finds a qualita- 
tively similar dependence of (C/ s )max(T c ) on the start- 
ing time To: The curves C max (T n ) move up as one in- 
creases To from 0.6 to l.Ofm/c (solid and dashed green 
lines). However, in contrast to the N-S initialization, the 
(C/ s )max(T c ) curves now show a strong dependence on 
relaxation time T n , rising monotonically with T n . The 
reason is that it takes some time for the bulk pressure 
II to develop large enough magnitudes to violate the 
positivity condition p+Tl > 0; again this happens typi- 
cally in regions where the matter is close to the phase 
transition. For larger relaxation times II moves away 
from its zero initial value more slowly, rendering the 
fluid more stable and resulting in a monotonic increase 
of (C/s)max(7c) with T n . For t h < 1 fm/c we find "univer- 
sal" ((/ s) max (T c ) — t h curves that do not depend on the 
shear viscosity rj/s (solid black, green and blue curves), 
but move upwards as we increase the starting time To. 
This is because the violation of the positivity condition 
p+II > then generally happens at early times t < 3 fm/c 
when the flow profiles are not yet significantly affected by 
shear viscous effects. For the two viscous fluid lines with 
r]/s = 0.08 and 0.16 (solid blue and green lines) one sees 
that they continue to overlap even for T n > 1 fm/c after 
they have broken away from the rj/s = line. In the 
ideal fluid (?y/s = 0) the phase transition generates large 
velocity gradients near the phase transition that generate 
locally large expansion rates, causing instability at lower 
values of £/s. Shear viscosity smoothes out these large 
gradients, as discussed in Ref. [20], allowing the fluid to 
evolve stably up to larger values of (/s. Bulk viscosity ( 
alone has no smoothing influence on sharp structures gen- 
erated by a phase transition. For a zero initial value for 
II, shear viscosity thus helps crucially in stabilizing the 
evolution of the viscous fluid against mechanical instabil- 
ities caused by strongly negative bulk viscous pressure, 
especially for large relaxation times T n . 

Very interesting is our finding that, for zero initial con- 
ditions, there is no limit on C max if one accounts for crit- 
ical slowing down of the bulk pressure dynamics near T c 
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via Eq. (1). In this case the bulk pressure, starting at 
zero, never grows sufficiently large to threaten mechani- 
cal stability of the fluid, irrespective of how large the bulk 
viscosity becomes at T c \ As the peak value {C/ S )(T C ) is 
increased, so is the time it takes II to evolve towards its 
Navier-Stokes value, and this never happens fast enough 
to violate the stability condition p+U > 0. 

Figure 6b shows the viscous entropy production for the 
maximally allowed bulk viscosities shown in Fig. 6a. Not 
surprisingly, viscous entropy production increases with 
shear viscosity rj/s and decreases when hydrodynamics 
is started later, with correspondingly smaller initial ex- 
pansion rates [21]. The dependence on (C/s) max (T c ) is 
non-monotonic, however. The reason is that the bulk 
viscous entropy production rate ~n 2 /(2£) depends not 
only on how large C is but also on how close II is to its 
Navier-Stokes limit, and this in turn depends on r n . 



VII. TOWARDS EXTRACTING rj/s FROM 
EXPERIMENTAL DATA: UNCERTAINTIES 
INTRODUCED BY BULK VISCOSITY 

Given the fact that bulk viscosity contributes to the 
viscous suppression of elliptic flow (see Fig. 5a), and as- 
suming that bulk and shear viscous effects cannot be sep- 
arated by studying other experimental observables, the 
question arises naturally how much of an irreducible un- 
certainty this will introduce into the extraction of the 
specific shear viscosity rj/s from experimental elliptic flow 
measurements. More precisely, if the QGP should turn 
out to be a "most perfect liquid" with "minimal" shear 
viscosity rj/s = with what kind of accuracy can we 

hope to verify this experimentally if bulk viscosity is the 
only quantity beyond our theoretical and experimental 
control? 

To answer this question, we used VISH2+1 to com- 
pute the differential elliptic flow of directly emit- 
ted pions (without resonance decay contributions) for 
200AGeV Au+Au collisions at 6 = 7fm, assuming the 
fireball medium to have constant specific shear viscosity 
?//s = l/47r but allowing the bulk viscosity £/s to vary 
over the entire range allowed by the mechanical stabil- 
ity criterium p+Tl > 0. In doing so we assumed a fixed 
shape of the temperature dependence of £/s as shown 
in Fig. 1 but let its normalization vary between C = 1 
and C max (r n ) where the latter is the largest value within 
the range of applicability of Israel-Stewart viscous fluid 
dynamics, shown in Fig. 6a. We allowed for two fixed 
values of 0.5 and 5fm/c for the bulk viscous relaxation 
time r n as well as for "critical slowing down" according 
to Eq. (1), and we studied both zero and Navier-Stokes 
initial values for the viscous pressure components. All 
calculations assume To = 0.6fm/c as starting time. The 
results are presented in Fig. 7 and Table I. 

Generically one observes that, even for minimal shear 
viscosity near the KSS bound, the shear viscous contri- 
bution to the elliptic flow suppression far exceeds the 




P T (GeV) 

FIG. 7: (Color online) V2(pr) for directly emitted pions from 
ideal and viscous hydrodynamics with Navier-Stokes (a) or 
zero (b) initial conditions for the viscous pressures. Shown 
are results for minimal shear viscosity t]/s = 1/Ati and bulk 
viscosities ranging from "minimal" (C = 1) to the maximal 
values from Fig. 6a that still allow for stable viscous evolution, 
for three choices of the bulk viscous relaxation time r„ . 



bulk viscous contribution. This is good news since it 
means that the uncertainty introduced into the extrac- 
tion of rj/s by theoretically poorly controlled bulk viscous 
effects remains limited and is, in fact, quite small, espe- 
cially if the real fireballs created in heavy-ion collisions 
do not completely saturate the KSS bound. On a more 
quantitative level, one finds that for pions with typical 
transverse momentum px — 0.5 GeV/c the elliptic flow is 
suppressed by just over 16% below the ideal fluid value if 
the expanding matter has only shear, but no bulk viscos- 
ity, and that this suppression increases to values between 
17% and 25% if bulk viscosity is added. The largest bulk 
viscous suppression is found for fixed relaxation times r n 
and zero initialization if the bulk viscosity is increased 
all the way up to its upper allowed limit. In these cases 
the additional suppression can be as large as 50% of the 
suppression found for the fluid with only (minimal) shear 
viscosity. If one takes into account that the evolution of 
the bulk viscous pressure slows down near T c where £/s 
is largest, the additional bulk viscous suppression never 
exceeds 20% of the shear viscous elliptic flow suppression, 
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TABLE I: Pion elliptic flow at pr = 0.5 GeV/c for b= 7fm 200 A GeV Au+Au collisions from ideal and viscous hydrodynamics, 
with different choices of initial conditions, bulk viscous relaxation times r n , and bulk viscosities (parametrized by C). The 
last of the 3 columns in each initialization block gives the viscous suppression of V2 at pr = 0.5 GeV/c in terms of the percent 
deviation from the ideal fluid baseline ( = 100%). 
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0.08 
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Eq. (1) 
Eq. (1) 


1 
100 


4.743 
4.656 


82.5 
80.9 


1 

1.3 


4.660 
4.615 


81.0 
80.2 



with 10-15% being a typical range (light blue and green 
curves in Fig. 7). 

An important caveat is, however, that for Navier- 
Stokes initial conditions the allowed maximal bulk vis- 
cosities are small, much below recent Lattice QCD esti- 
mates [76]. If larger values are realized by Nature, they 
invalidate the use of viscous hydrodynamics, at least at 
early times [89-92] . The problems in this case arise from 
the large bulk viscosity in a thin layer near the transverse 
edge of the fireball where the matter is close to T c . It is 
only in this region that the viscous hydrodynamic de- 
scription breaks down. Since the problematic factor, the 
scalar expansion rate 9, decreases initially very rapidly, 
these initially unstable fluid regions move quickly back 
to mechancal stability. Since the momentum anisotropy 
does not develop instantaneously, we find it hard to be- 
lieve that the existence of this unstable external layer has 
much influence on the evolution and final value of the el- 
liptic flow, and one should get very similar results from 
simulations in which the initial bulk viscous pressure II 
is restricted by hand to values below the threshold for vi- 
olating the positivity condition p+H > 0. If this is indeed 
the case, the results presented in this Section show that 
bulk viscosity, even if theoretically not well controlled, 
will not introduce large uncertainties into the extraction 
of r)/s from elliptic flow data. 

VIII. CONCLUDING REMARKS 

The present study shows that bulk viscosity, as long 
as it is small enough that in expanding heavy-ion colli- 
sion fireballs the negative bulk viscous pressure does not 
become larger than the thermodynamic pressure, affects 
the elliptic flow of the final hadrons much more weakly 
than does shear viscosity. So, as long as the expanding 
fireball can be described by viscous fluid dynamics, it is 
possible to extract its shear viscosity (even if it is as small 
as -j| K gg = j^) with good accuracy from a comparison of 
viscous hydrodynamic simulations with experimental el- 



liptic flow data. Accounting for the critical slowing down 
of viscous bulk pressure dynamics near T c , we showed 
that any contamination from bulk viscosity Q/s is < 20% 
(for much of the parameter space it is even < 10%), and 
that its relative importance decreases further if rj/s is 
larger than the KSS bound. 

However, we also saw that the stability condition 
p+H > is very restrictive and easily violated if the peak 
value of Q/s near T c reaches values close to those esti- 
mated from Lattice QCD [76] and from some strong cou- 
pling approaches [66], and if the bulk viscous pressure 
IT approaches its Navier-Stokes limit 11 = — Qd-u. When 
this occurs (typically at early times when the scalar ex- 
pansion rate is largest, in a thin layer around T c close 
to the transverse edge of the fireball), the viscous fluid 
dynamical description breaks down. Our analysis shows 
that the phenomenon of "critical slowing down" can play 
a crucial role in preventing this from happening. Kinetic 
theory for weakly coupled systems [56, 93] and a recent 
analysis by Buchel of strongly coupled systems [66] sug- 
gest that the same microscopic physics (namely growing 
correlation lengths due to critical fluctuations) that gen- 
erate a peak of Q/s at T c also causes the relaxation time 
T n for the bulk viscous pressure to grow and possibly di- 
verge at T c even while Q/s itself remains finite. When 
using the model Eq. (1) for a temperature dependent r n 
inspired by these ideas we saw that, unless II is initialized 
at its Navier-Stokes limit, it never reaches it during the 
short time span of a heavy-ion collision in those fireball 
region where Q/s peaks and II could thus become very 
large. This reduces the problem of applicability of vis- 
cous hydrodynamics at early times to a question of initial 
conditions for II, especially in that thin transverse layer 
where (after local equilibrium is reached) the tempera- 
ture happens to be close to T c . 

Determining these initial conditions (as opposed to 
guessing them as we have done here) requires a theoreti- 
cal description of the early pre-equilibrium evolution and 
Landau-matching the corresponding energy-momentum 
tensor to its viscous fluid dynamic form, Eq. (2) (in the 
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spirit of Ref. [92] but generalized from 0+1 to 2+1 di- 
mensions). At this point we lack the tools for doing this. 
Let us, however, make a few comments in anticipation of 
completion of that task. Consider a small fireball region 
that is just reaching local thermal equilibrium at a tem- 
perature close to T c and undergoing self-similar boost- 
invariant longitudinal expansion while transverse expan- 
sion is negligible. Let us also assume that at this point 
in time the bulk viscous pressure in the region is large 
and negative, leading to negative effective total isotropic 
pressure and causing the fluid to be mechanically unsta- 
ble. What will happen? The fluid will begin to rupture, 
forming little voids, and if the region were to remain in a 
state of negative total pressure, it would eventually frag- 
ment. However, since the considered region is undergoing 
rapid expansion and cooling, it will quickly exit from its 
state of mechanical instability. Furthermore, during the 
short period of instability the hydrodynamic growth of 
voids will be hampered by the large value of the relax- 
ation time T n . By the time the considered region becomes 
mechanically stable again, we expect it to be riddled with 
small holes, but otherwise intact. The small voids formed 
during the period of instability will re-collapse by cavita- 
tion, and the region will quickly re-equilibrate due to the 
now much shorter relaxation time below T c . No whole- 
sale breakup of the fluid will occur, due to lack of time. 
Similar arguments hold later when the bulk of the matter 
in the center of the fireball passes through T c , only that 
in this case the viscous bulk pressure may never grow 
large enough to generate mechanical instability, due to 



critical slowing down. 

In summary, unlike the authors of Ref. [48], we do 
not expect any dramatic macroscopic phenomena trig- 
gered by the transient mechanical instability arising from 
possibly large, but short-lived negative bulk pressures in 
fireball regions passing through the hadronization phase 
transition. For this reason we believe that a modified vis- 
cous hydrodynamic treatment, where one limits by hand 
the growth of the viscous bulk pressure so that it al- 
ways remains below the instability threshold [94], will 
not lead to impermissible distortions of the real (non- 
equilibrium) dynamics in the (small) space-time regions 
whose description lies outside the hydrodynamic domain. 
This is important for future viscous hydrodynamic stud- 
ies of heavy-ion collisions with fluctuating and granular 
initial conditions [95] which are more realistic than the 
smooth initial profiles presently used. 
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